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ABSTRACT 

Accurate neutrino transport has been built into spherically symmetric simulations of stellar core col- 
lapse and postbounce evolution. The results of such simulations agree that spherically symmetric 
models with standard microphysical input fail to explode by the delayed, neutrino-driven mecha- 
nism. Independent groups implemented fundamentally different numerical methods to tackle the 
Boltzmann neutrino transport equation. Here we present a direct and detailed comparison of such 
neutrino radiation-hydrodynamical simulations for two codes, AGILE-boltztran of the Oak Ridge- 
Basel group and vertex of the Garching group. The former solves the Boltzmann equation directly 
by an implicit, general relativistic discrete angle method on the adaptive grid of a conservative implicit 
hydrodynamics code with second-order TVD advection. In contrast, the latter couples a variable Ed- 
dington factor technique with an explicit, moving-grid, conservative high-order Riemann solver with 
important relativistic effects treated by an effective gravitational potential. The presented study is 
meant to test both neutrino radiation-hydrodynamics implementations and to provide a data basis for 
comparisons and verifications of supernova codes to be developed in the future. Results are discussed 
for simulations of the core collapse and post-bounce evolution of a 13 Mq star with Newtonian gravity 
and a 15 M star with relativistic gravity. 

Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics — relativity — 
methods: numerical 



1. INTRODUCTION 

Computer simulations are becoming more and more 
popular. They allow investigations of physics on office 
desks rather than explorations through hands-on experi- 
ments (This does not imply a transition from hard work 
to gaming). In the industrial context, the two approaches 
are not separable: the computer codes have to be val- 
idated. After a computer design has been completed, 
its relation to reality will inevitably be assessed in the 
manufacture and evaluation of prototypes. How about 
the growing importance of computer simulations in as- 
trophysics - where are the measurements found to test 
aspects of a complex computer code in idealized setups, 
and where are the prototypes that validate the quality 
of the results in the targeted application? The first step 
of code development is accompanied by the verification 
of partial aspects of the code in simplified test problems 
where the solution is analytically known. The code can 
be improved step by step. Additionally, laboratory ex- 
periments may be used to further verify the code with 
accurate measurements in idealized setups. The transi- 
tion to code validation is made when its capability to 
handle complicated coupled processes is tested and the 
completeness of the physical description is evaluated. In 
the industrial context, this is achieved by more compre- 
hensive experiments and measurements, or, ultimately, 
by the comparison of computer designs with the prop- 

1 CITA, University of Toronto, Toronto, Ontario M5S 3H8, 
Canada 

2 Physics Division, Oak Ridge National Laboratory, Oak Ridge, 
Tennessee 37831-6354 

3 Department of Physics and Astronomy, University of Ten- 
nessee, Knoxville, Tennessee 37996-1200 

4 Max-Planck-Institut fiir Astrophysik, Karl-Schwarzschild- 
Strasse 1, D-85741 Garching, Germany 



erties of manufactured prototypes. We would now be 
tempted to relate the validation of computer generated 
results with real life prototypes in manufacturing to the 
comparison of astrophysical simulations with astronom- 
ical observations. This would, however, circumvent the 
goal of astrophysical simulations: One does not assume 
unknown physics in industrial design. Perfect agreement 
between a computer simulation and the behavior of a 
prototype can indeed be seen as proof of the quality of 
the computer code. The situation is different in astro- 
physics, where the understanding of the physics of an 
event is rather the goal than the ingredient. The com- 
parison between simulation and observation is essential 
to demonstrate the physical understanding of the event, 
it cannot at the same time be used to qualify code per- 
formance. The gap in code validation between detached 
analytical test calculations and the astroph ysical applica- 
tion c an be bridged by code comparisons l| C alder et al. I 
2002), based on the assumption that different numeri- 
cal approaches are likely to show different strengths and 
weaknesses in simulations of complex physical systems. 
Differences in the simulation results are an indicator for 
uncertainties in the numerical methods. 

In the present paper we document the detailed com- 
parison of results from different supernova codes. Both 
of our codes aim to provide a solution to the Boltzmann 
neutrino transport equation in spherical symmetry. 
This is achieved by fundamentally different numerical 
methods: The code of the Oak Ridge-Basel group 
(agile-BOLTZTRAN) consists of a general relativistic 
time-implicit discrete-angle (Sat) Boltzmann solver, 
which is coupled in an operator split fashion to a general 
relativistic time-implicit hydrodynamics solver with a 
dynamical adaptive grid. It implements a direct finite 
difference representation of the Boltzmann equation 
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iMezzacappa fe Bruenn I Il993al iMezzacappa fc Messer I 
1999t ILiebendorfer. Rosswog. & Thielemann I l2002t 
Liebendorfer et al. Il2004j) . The Garching code, vertex , 
is a one-dimensional version of a program that was 
developed to perform multi-dimensional supernova 
simulations with accurate ray- by-ray neutrino transport. 
It is based on the explicit, moving-grid, finite-volume 
hydrodynamics code PROMETHEUS, which employs a 
Riemann solver for constructing the solution of the 
hydrodynamics equations. The neutrino transport is 
handled in an operator-split step and is calculated by 
a variable Eddington factor closure of neutrino energy, 
number, and momentum equations, where the variable 
Eddington factor is derived from the formal solution 
of a spherically averag ed model Boltzmann equation 
ijRampp fe Jankal [20021. In spherical symmetry, there 
is only one "ray" for the solution of the moments 
equations and no spherical averaging is necessary for the 
model Boltzmann equation. Therefore, a convergence 
in the iterations between the moments equations and 
the closure from the model Boltzmann equation pro- 
vides in spherical symmetry a solution of the complete 
Boltzmann equation. 

This work has two main goals, (i) The direct compar- 
ison of two codes applied to the same challenging astro- 
physical scenario with concerted physics (spherical sym- 
metry, progenitor models, nuclear and weak interaction 
physics, general relativistic effects), (ii) The production 
of reference results to test future supernova codes in the 
spherical limit. Machine-readable data files are included 
in the electronic edition of the Journal. 

Neutrinos play a crucial role in collapsing cores of mas- 
sive stars. The loss of electron lepton number by the 
production and escape of electron neutrinos determines 
the collapse dynamics and the position where the su- 
pernova shock forms. Energy and lepton number trans- 
port by neutrino diffusion also govern the evolution of 
the nascent neutron star. The energy transfer by neu- 
trinos to the medium that surrounds the protoneutron 
star may revive the sta lled accretion front and thus drive 
a de layed explosion jWilson ( |1985t IBethe fe Wilson 1 
1985). * Neutrino interactions moreover set the proton- 
to-nucleon ratio and therefore the conditions for nu- 
cleosynthesis in the innermost supernova ejecta. Sus- 
tained energy deposition near the protoneutron star sur- 
face causes a post-explosion outflow of baryonic matter, 
the so-called neutrino-driven wind, which is discussed 
as a p otential site for the formation of r-process ele- 
ments (IWooslev et al. ll994tlTakahashi. Witti. fc Janka I 
[1991 ISumivoshi et al. I 120001 IWanaio et al. I 120011: 
iThompson. Burrows, fc Mever 112001^1 . 

Deep inside the protoneutron star the absorption and 
scattering mean free paths of neutrinos are very small 
and therefore neutrinos diffuse and are in chemical equi- 
librium with the stellar plasma. With decreasing density 
the neutrino interaction lengths become larger, before, fi- 
nally the neutrinos can stream freely. Since the reaction 
cross sections rise steeply with the neutrino energy, low- 
energy neutrinos decouple from the stellar background at 
higher densities. Most electron flavor neutrinos emerge 
from the accreting material at the base of the cooling re- 
gion under semi-transparent conditions and propagate to 
the heating region where their angular distribution influ- 
ences the energy deposition behind the accretion front. 



Neither diffusion nor free streaming is a good approxi- 
mation in this important region where neutrinos strongly 
couple the dynamics of different layers on a short prop- 
agation time scale. 

An accurate treatment of the neutrino transport 
and of neutrino-matter interactions therefore requires 
the combination of neutrino sources at one loca- 
tion with neutrino opacities at other locations as de- 
scribed by the energy- and angle-dependent Boltz- 
mann transport equation. The solution of the Boltz- 
mann equation is also desirable to test approxima- 
tions, the most elaborate of which are certainly multi- 
group flux-limited diffusion llBru enn I IT9 851 iMvra et al I 
Il987t iBruenn. DeNisco. fc Mez zacappa II2001D and two- 
mom e nt closure schemes l|B ludman fc Cer nohorskvl 
119951 ISmit. Cernohorskv. fc Dullemond 1 11997ft . With 
the growing computer capability it has become feasible 
to provide solutions to the Boltz mann transport equation 
not on ly for the collapse phase ({Mezzacappa fc Bruenn I 
1993a), but also in consistent dynamical simulations 
of the post-bounce evolu ti on llRamp p fc Ja nka 1 120001 
Mezzacappa et al.~l I200H ILiebendorfer et al. I 120011: 
Thompson. Burrows, fc Pinto II2003D . 

The paper is organized as follows. We will describe 
in Sect. |2 the stellar models and the physical ingredients 
that constitute the problem to be solved by our numerical 
methods. In Scct.|21we will briefly resume characteris- 
tic features and capabilities of both neutrino radiation- 
hydrodynamics codes. In Sect. 0] our results for the two 
considered stellar models will be discussed with special 
focus on the differences between the runs. In Sect. [S] we 
shall summarize our findings and draw conclusions. 

2. DESCRIPTION OF THE MODEL 

2.1. Progenitors 

We present in this paper two different models. One 
of them includes only the most essential physical ingre- 
dients; e.g., the transport of electron flavor neutrinos 
and antineutrinos, but not the heavy-lepton neutrinos 
and antineutrinos. It is meant to serve as a guideline 
for future code development and transport approxima- 
tions. The model is based o n a 13 M Q progenitor of 
iNomoto fc Hashimoto 1 l|1988fl . The 13 M progenitor 
model has a tradition in supernova simulations, its ex- 
ceptionally small iron core sustained the hope to pro- 
duce prompt explosions. We call this run N13, as it is 
based on Newtonian gravity. A second run was launche d 
from a 15 M Q progenitor of IWooslev fc Weaver I |l995). 
This progenitor has been widely used in supernova in- 
vestigations, as it provides a model of a massive star in 
the middle of the mass range that is expected to end its 
life in a supernova. The run takes into account all neu- 
trino flavors w ith "standard" input physics as listed in 
IBruenn I (£l985). General relativistic effects are included 
in this physically more complete run, G15, and the input 
physics has been extended by ion-ion correlations and 
nucleon-nucleon bremsstrahlung . 

2.2. Radiation Hydrodynamics in Spherical Symmetry 

The stellar progenitor model is evolved in time by 
means of the hydrodynamics and neutrino transport 
equations. Since there is no danger of grid entan- 
gling in spherical symmetry, we make use of the free- 
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dom in AGILE-BOLTZTRAN to choose orthogonal comov- 
ing space-time coordinates and describe the interactions 
in the collision term in the most convenient comov- 
ing frame for the neutrino f our-momentum (see e.g., 
iCardall fc Mezzacaooa 1 ((2003) for a generalized discus- 
sion of coordin ate choices for the radia t ion transport). 
In the metric of iMisner fc Sharp I (| 1 964(1 : lMav fc White I 
(196§, 



da 2 +r 2 (dtf 2 + sin 2 ddtp 2 ) 



ds 2 = —a 2 (cdty 



1 dr s 
fda, 

the equations of hydrodynamics in spherical 
symmetry in the presence of a radiation field 
the 



be written in the form jLindauistl [l966; 
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The metric is based on an enclosed baryon mass label, 
a, and a coordinate time, t. It refers to the areal radius, 
r, the "Lorentz" factor, T = yl-j- (u/c) 2 — 2Gm/ (rc 2 ), 
and the lapse function, a. The angles $ and tp describe a 
two-sphere. The quantity u = or x drjdt takes the place 
of a fluid velocity and m is the enclosed gravitational 
mass, proportional to the enclosed total energy. Both of 
our codes split these equations into hydrodynamics equa- 
tions and transport equations, vertex employs Eulerian 
coordinates which can be obtained by a coordinate trans- 
formation. The fluid is specified in its rest frame by the 
rest mass density, p, specific internal energy, e, and elec- 
tron number fraction, Y e , for which an additional evo- 
lution equation (lepton number conservation) is solved. 



The hydrodynamics equations are closed by the equa- 
tion of state, which gives the pressure, p, as a function 
of density, internal energy, and electron number. The 
zeroth, J, first, H j 'c, and second, K, angular moment of 
the monochromatic neutrino intensity (normalized by the 
rest-mass density) are determined by the transport equa- 
tion which includes the interactions listed in Table Urn 
the collision integral. Eq. (0) determines the lapse func- 
tion. Eq. jnjl allows us to integrate outward from a = 
to obtain the total energy; it translates to the Poisson 
equation in the Newtonian limit. Eq. (jSJ) defines the gen- 
eral relativistic analogue to the Newtonian volume. Eq. 
(@J describes the change of radial momentum; to leading 
order it is proportional to —Gm/(cr) 2 + (J — 3K)/r. Eq. 
(0 describes the evolution of the total energy. Note that 
there are no contributions from terms of zeroth and first 
order in c. Eq. J2J) relates the evolution of the specific 
volume to the divergence of the velocity field. 

2.3. Nuclear and Weak Interaction Physics Input 

The equation of state describes the thermodynamic 
state of a fluid element based on density, p, tempera- 
ture, T, and the composition. The relation between the 
specific internal energy and the pressure closes the sys- 
tem of hydrodynamics equat ions. For this comparison w e 
use the equation of state of lLattimer fc Swestvl <|1991|) . 
It assumes nuclear statistical equilibrium and is based 
on a liquid drop model for a representative nucleus with 
atomic number A and charge Z, surrounded by free al- 
pha particles, protons, and neutrons. These baryons are 
immersed in an electron and positron gas that equili- 
brates with a photon gas by the pair creation process. 
The incompressibility modulus can be adjusted. We use 
a value of K — 180 MeV. Above nuclear density, where 
no isolated individual nuclei are present, the transition 
to a proton-neutron-electron gas is made by a Maxwell 
phase transition. In any of these cases, the nuclear com- 
position at given temperature and density is determined 
by the specification of the electron fraction Y e , At low 
densities/temperatures the vertex code switches to an 
equation of state that considers electrons, positrons, pho- 
tons, nucleons and nuclei consis tent with the compositio n 
used in the progenitor model l|Rampp fc Janka ^ 2002). 
The switch is triggered in N13 and G15 by a density 
threshold of p < 6 x 10 7 g/cm 3 . AGILE-boltztran ap- 
plies the same switch in the run N13, but considers in the 
low density domain only one nucleus, 28 Si. In run G15, 
silicon is converted to NSE under en ergy conservation at 
a bur ning temperature of 0.44 MeV ( Mczzac appa et al. I 
I2001 . 

As for the neutrino-matter interactions, we have cho- 
sen to use the conventional ("standard") opacities, i.e., 
a descrip tion w hich follows closely the o ne detailed by 
iBruennl l(T985|) : iMezzacappa fc Bruenn I l|1993c|) . Note, 
however, that there are small differences in the neutrino 
description employed by the two groups. While AGILE- 
BOLTZTRAN treats the p/r neutrinos and antineutrinos 
separately, they are combined to one species in vertex. 
In order to save computer time, usually only electron 
neutrinos are considered in the vertex calculations dur- 
ing the core collapse phase. Tests have shown that taking 
into account also electron antineutrinos and the heavy- 
lepton neutrinos leads to only minuscule differences in 
this phase of the supernova evolution (their inclusion in 
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Table 1 

Overview of all neutrino-matter interactions 
considered in the n13 runs. 



Reaction™ Reference 



uA 




vA 


Bruenn (1985) (no ion-ion correlations!) 








i/N 


Bruenn (4985); Mczzacappa & Bruenn 


1993c) 


i/ e n 




e - p 


Bruenn (1985); Mczzacappa & Bruenn 


1993c) 


i'eP 




e+n 


Bruenn (1985); Mczzacappa & Bruenn 


1993c) 


u e A' 




e~A 


Bruenn (1985): MezzacaDDa &: Bruenn 


1993c) 


vv 




c~c+ 


Bruenn (1985) 










Mczzacappa & Bruenn (1993c); 





^crnohorsky (1994) 

a In the first column the symbol u represents any of the neutri- 
nos v c , P c (heavy-lepton neutrinos are neglected in N13). The 
symbols e _ , e + , n, p and A denote electrons, positrons, free neu- 
trons and protons, and heavy nuclei, respectively, the symbol N 
means n or p. The references point to papers where informa- 
tion can be found about the approximations employed in the rate 
calculations. Details about the numerical implementation can be 
found in the methodical pap e rs oflR ampp &: Janka l 120021) and 
Mczzacappa & Mcsser (1999); Licbcndorfcr et al. (2004), respec- 
tively. The G15 runs additionally include reactions with fi- and 
t- neutrinos and their antiparticles, ion-ion correlation effects in 
neutrino-nuclei interactions, and nucleon-nucleon bremsstrahlung. 

the postbounce phase, however, is important as demon- 
strated by the comparison between models N13 and G15 
in Sect. 0}. The Garching group routinely includes 
nucleon-nucleon bremsstrahlung as a source (or sink) for 
neutrino-antineutrino pairs. The decrease of coherent 
neutrino scattering off heavy io ns due to ion-ion cor- 
relations and electron screening JItoh 1119751 iHorowitz I 
Il997t iBruenn fc Mezzacappa Il997[) is also taken into ac- 
count. These improvements have been switched off in 
the N13 models and added to agile-boltztran such 
that they are consistently included in both G15 runs. 
In both codes, the implementation of the ion-ion correla- 
ti qns has been updat ed with the structure function given 
in lltoh et al. 1 1)2004|) . Because of the large mass contrast 
between species with A < 4 and the representative heavy 
nucleus we omit the somewhat arbitrary averaging of the 
effect over species and consider only the representative 
nucleus for the calculation of the ion separation. 

3. NUMERICAL METHODS 

The two codes follow very different approaches to eval- 
uate the radiation moments J, H, and K. Contrary to 
flux-limiting and "gray" transport methods, neither of 
our methods needs to make assumptions about the an- 
gular or the spectral distribution of the radiation field. 
Important features and implementation details of the two 
codes are summarized in the following subsections. 

3.1. AGILE-BOLTZTRAN 

The concept and first imp lementation of BOLTZ- 
tran ha s been developed in iMezzacappa fc Bruenn I 
l)1993albllcT) for core collapse simulations in the order 
v/c approximation. Essential for the computational 
efficiency of the implicit scheme is the storage of the in- 
teractions in a dynamical table which delivers consistent 
derivatives of all cross sections and thermodynamical 
quantities for the Newton- Raphs on iterations in the solu- 
tion of the nonlinea r equations ( Mczzac appa fc Messer I 
I1999L iMesser l l2000). For the highly dynamical situation 
after bounce, BOLTZTRAN has been coupled to the 



hydrodynamics code agile. The finite differencing has 
been revised for energy conservation an d extended to 
solve th e general relativistic Eas. (12151 (ILiebendorfer I 
20001 I Liebendorfer R osswog. fc Thielemann I 120021 
Liehendorfer et al. II2004D. 

3.1.1. Hydrodynamics 

The hydrodynamics part of the Lagrangian Eqs. (121711 
is solved by implicit conservative finite differencing. One 
strength is the implement ation of a dynamically moving 
adapt ive g rid following | Winkler. Norman, fc Mihalas I 
1)19841) and iDorfi fc Drurv I l)1987|) . In the general rel- 
ativistic extension, it is equivalent to a resolution- 
dependent choice of shift vectors, that allow a continuous 
coordinate translation in the radial direction during the 
evolution of the model. Artificial viscosity has been in- 
cluded in a consistent, but ca usality violating manner 
based on the tensor viscosity of lTscharnuter fc Winklerl 
(1979). It provides the mechanism for energy dissipation 
in the shock front and defines the shock width such that 
the number of attracted adaptive grid points does not 
grow beyond limits. A major advantage of the adaptive 
grid is the dynamical allocation of computational zones 
to regions where they are needed. The zoning for hy- 
drodynamics and neutrino transport is always kept con- 
gruent for consistency reasons. During the evolution of 
the model, one group of grid points follows the accretion 
front, while another resolves the steep density gradient 
between the outer layers of the protoneutron star and the 
infalling matter, where most of the electron flavor neu- 
trinos stream away. The use of artificial viscosity is not 
a disadvantage of the method, as it only plays a minor 
local role in stabilizing the shock front (the shock width 
is set to a few percent of the shock radius and captured 
by ~ 10 moving grid points). A disadvantage of the 
adaptive grid approach in an earlier imp l emen tation (see 
ILiebendorfer. Rosswog. fc Thielemann (2002) for a de- 
tailed description) is the numerical diffusion introduced 
by the first-order advection in regions where the adaptive 
grid does not concentrate its zones, e.g. in a rarefaction 
wave or at sharp discontinuities of the composition. Dur- 
ing this comparison, improvement has been achieved by 
upgrading agile with second order total variation di- 
minishing (TVD) advection based on a Van Leer flux 
limiter. 

3.1.2. Neutrino Transport 

The neutrino transport part, boltztran, solves the 
Boltzmann transport equation, 

acot a oa 
\r a dr) dfj, ^ ^ ^ ^ 

V aedt rc ) d P 

in a finite difference representation implementing the dis- 
crete ordinates, or Sat, method. The evolved quantity 
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is the neutrino distribution function, F(t, a, /i, E), as a 
function of time t, rest mass a enclosed in a sphere at ra- 
dius r, propagation angle cosine [i with respect to radial 
direction, and neutrino energy E. Neutrinos in specific 
beams are created and destroyed according to the col- 
lision term, C, which includes the interactions listed in 
Tabled It is assumed that the neutrinos propagate freely 
between interactions. The free-particle motion along 
geodesies between collisions introduces the many correc- 
tion terms apparent on the right hand side of Eq. l|5)l. 
They stem from the use of spherical coordinates in com- 
bination with a description of the neutrino phase space in 
a comoving frame. Nevertheless, all terms can be labeled 
with a physical effect. In order of appearance in Eq. JSJ 
these are, the time derivative of the neutrino distribution 
function, the propagation of neutrinos, the angle correc- 
tion due to neutrino propagation, the angular aberration 
correction due to observer motion, the frequency shift in 
the gravitational potential, and the Doppler frequency 
shift due to observer motion. It is essential for the suc- 
cessful finite difference representation of Eq. iJSJ that 
it is upward compatible with simple special cases of the 
transport equation. Basic physical properties can be de- 
termined by the evaluation of expectation values with 
the neutrino distribution functions for various operators. 
The finite difference representation should support, e.g., 
the diffusion limit, total energy conservation, and the 
conservation of lepton number. 

3.1.3. Parameter settings 

Both runs, N13 and G15, were performed with 103 
adaptive spatial zones ranging from the center of the 
progenitor star to about 7000 km. A constant pressure 
boundary condition was used at the barely moving sur- 
face. The neutrino energy was resolved with 20 geomet- 
rically increasing energy groups, the first centered at 3 
MeV and the last at 300 MeV. The propagation angle 
has been discretized with 6 angles suitable for Gaussian 
quadrature. Roughly 3000 time steps have been used for 
collapse and 7000 for the postbounce phase. The run 
N13 has been evolved with an order v/c approximation 
of Eqs. 1)218(1 and run G15 with the general relativistic 
equations. 

3.2. VERTEX 

Independently from the efforts of the Oak Ridge- 
Basel collaboration the Garching supernova group has 
treated the Boltzmann transport problem for neutri- 
nos in core-collapse sup ernovae with a new variable Ed- 
dington fac tor method ( Rampp 200(J IRampp fc Jankal 
120001 12002ft . and has coupled it to the PROMETHEUS 
hydrodynamics code. T he combined program allows 
for spherically symmetric l|Rampp fc Janka J [2000jj200^l 
as well as multi-dimensio nal simulations ijJanka et al. I 
12004 iBuras et al. 1 1200!) . The spherically symmetric 
"core" of the program, which was used for the calcu- 
lations described below, will be referred to by the name 
vertex (Variable Eddington factor Radiative Transfer 
for supernova EXplosions). 

3.2.1. Hydrodynamics 

The integration of the equations of hydrodynamics 
is performed with the Newtonian finite-volume code 



PROMETHEUS l|Frvxell. Miiller. fc Arnett I I1989J) which 
was s upplemented by additional problem specific fea- 
tures (|Keil Ill997|) . PROMETHEUS is a direct Eulerian, 
time-explicit imple mentation of the Piecewise P arabolic 
Method (PPM) of lOolella fc Woodward I lfl98l . As a 
Godunov scheme of third order in space and second-order 
in time with a Riemann solver, it is particularly well 
suited for following discontinuities in the fluid flow like 
shocks or boundaries between layers of different chemi- 
cal composition. A notable advantage is its capability of 
tackling multi-dimensional problems with high computa- 
tional efficiency and numerical accuracy. The code makes 
use of t he "Consistent Multiflu id Advection (CMA)" 
method l)Plewa fc Miiller 1 Q.999) for ensuring accurate 
advection of different chemical components in the fluid, 
and switches from the o riginal PPM me thod to the more 
diffusive HLLE solver o flEinfeldt I l)1988|) in the vicinity of 
strong shocks to avoid spurious oscillations (the so-called 
"odd-even decoupling" phenomenon) when such shocks 
are aligned with one of the coordinate lines in multi- 
dimensional simulation s l)Quirk 1 119941: [Kifonidis J l2000t 
IPlewa fc Mime71l200l . 

3.2.2. Neutrino transport 

The variable Eddington factor scheme for the neutrino 
transport, its coupling to the hydrodynamics part, and 
application to a numb er of test problems is de scribed in 
much detail elsewhere l|Rampp fc Janka"l f2002). Here we 
will only briefly summarize the characteristic features of 
the method. The couple d set of equations of hyd rody- 
namics (Eqs. (l)-(4) in IRampp fc Janka I p002n and 
radiation transport (Eqs. (6)— (8) ibidem) is equivalent 
to Eqs. ©-0 m the order (v/c) limit. The equa- 
tions are also split into a "hydrodynamics part" and a 
"neutrino part" and solved independently in subsequent 
("fractional") steps. But the hydrodynamics and the 
transport solver can use radial grids and/or time steps 
that are different from each other. 

In the neutrino transport method the integro- 
differential character of the Boltzmann equation is tamed 
by applying a variable Eddington factor closure to the 
neutrino energy and momentum equations (and the si- 
multaneously integrated first and second order moments 
equations for neutrino number). For this purpose the 
variable Eddington factor is determined from the formal 
solution of the Boltzmann equation on so-called "tan- 
gent rays". They coincide with radiation characteristics 
in Newtonian geometry. The system of the Boltzmann 
equation and its moments equations is iterated until con- 
vergence is achieved. The integration of the transport 
equations is implicit in time. 

General relativis tic effects are treated only approxi- 
mately in the code l|Rampp fc Janka l2002j) . The current 
version contains a modification of the gravitational po- 
tential by including correction terms due to pressure and 
energy of the stellar medium and neutrinos, which are 
deduced from a comparison of the Newtonian and rela- 
tivistic equations of motion. The neutrino transport con- 
tains gravitational redshift and time dilation, but ignores 
the distinction between coordinate radius and proper ra- 
dius. This simplification is necessary for coupling the 
transport code to the basically Newtonian hydrodynam- 
ics. We shall demonstrate in this paper that these ap- 
proximations work satisfactorily well for the core collapse 
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and the early postbounce phase (see Sect. 14.3(1 . Moder- 
ate quantitative but no qualitative differences from the 
full relativistic treatment are mainly found at late times 
after the accretion front has started to retreat. 

3.2.3. Parameter settings 

For the neutrino transport in the N13 run an Eule- 
rian radial grid with 235 radial zones (255 tangent rays) 
spaced logarithmically between and 10 000 km was 
used. The neutrino spectrum between and 380 MeV 
was discretized with 21 geometrically zoned energy bins, 
the center of the first bin being located at 2 MeV. The 
hydrodynamics, on the other hand, was solved on a radial 
grid of 400 zones which are moved with the stellar fluid 
during core collapse. Shortly after core bounce both ra- 
dial grids were rezoned such that inside of a radius of 400 
km the zones of the transport grid coincide with those of 
the hydro grid. For the post-bounce evolution the coor- 
dinates of the hydrodynamics grid (as well as those of the 
transport grid) remained fixed in time. Concerning the 
resolution and definition of numerical grids the same pa- 
rameters were chosen in the G15 run with the exception 
that 19 energy bins between and 380 MeV were used 
instead of 21 groups and a rezoning of the radial grid 
was necessary at ~ 200 ms after bounce. The new grid 
(hydrodynamics and neutrino transport) employs 40 ad- 
ditional radial zones in order to adequately represent the 
steepening density gradients at the surface of the nascent 
neutron star. 

4. COMPARISON OF THE RESULTS 

The simulations produce a sizeable amount of data, 
even if they are confined to one spatial dimension. Hence, 
we start with an overview of what we are going to com- 
pare and how the comparable quantities are derived from 
the code-specific results. During collapse and bounce it is 
quite natural to choose the enclosed baryon mass as spa- 
tial coordinate. It labels individual fluid elements and al- 
lows one to trace the history of each fluid element. Some 
tenths of milliseconds after bounce, the neutron star be- 
comes rather static with mass accretion being essentially 
stationary. Therefore the presentation of the quantities 
as functions of radius is more appropriate. We present 
the three independent quantities that determine the ther- 
modynamic state of the fluid element as measured by an 
observer comoving with the fluid: the rest mass den- 
sity, the entropy per baryon, and the electron fraction. 
Furthermore, the radial velocity of the fluid element is 
displayed. 

Note that in case of the relativistic model, G15, the co- 
ordinate independent change in areal radius per proper 
time, u = a~ 1 dr/dt 1 is plotted in the (b)-panels of 
Figs. OD El El and El for the AGILE-BOLTZTRAN run, 
whereas the velocity, v — dr/dt, is unchanged from the 
Newtonian case in the general relativistic approximation 
of vertex. Typical deviations of the metric coefficients 
from unity are shown in Table In summary, these de- 
viations are of order 1% in the preshock region, increase 
from 1% to 6% at the neutrinosphere during the sim- 
ulation, and reach maximum values around 20% at the 
center of the star in case of the lapse function and around 
10% close to the center in the case of T. 

The neutrino transport quantities are represented by 
the neutrino luminosity and rms energy profiles as mea- 
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"Listed are the lapse function, a, and the metric coefficient, T, 
by their deviations from unity in percent at several post-bounce 
times tpj, during the G15 simulation with AGILE-BOLTZTRAN. The 
indices of a and T refer to chosen locations where the values are 
given. The index s corresponds to the position outside of the shock 
or accretion front, the index v corresponds to the neutrinosphere 
of the heavy lepton neutrinos, and the index max stands for the 
maximum value in the whole star. 



sured in the fluid rest frame. We also discuss selected 
quantities as functions of time, e.g. the trajectories of 
fluid elements, the position of the (accretion-)shock, the 
conditions at the center, or the neutrino luminosities and 
rms energies sampled at 500 km radius (as an approxima- 
tion to radial infinity). The physical time in both runs 
is synchronized at bounce (t — 0), the moment when 
the central density reaches a local maximum immedi- 
ately before the shock is formed. Negative times in the 
simulations therefore point to instances before bounce. 
We define the shock position, R s , by the maximum of 
the velocity divergence (i.e. maximum compression, cf. 
Eq. Q), 



x(R s ) — max [x(r)] , x = 



d(4-Kr 2 u) 
dV ' 



The luminosities and rms energies are given in the co- 
moving frame as a function of the neutrino phase-space 
distribution function F(t, a, [i, E) according to 



L(t, a) = Anr 2 p 



2ttc 
Jhcf 



J F(t,a,n,E)E 3 dEndfi 



(e(t,a)) r 



' / F(t,a,n,E)E' i dEdp 
J F(t, a, n, E)E 2 dEdp ' 



In the general relativistic case, the definition of 
F(t,a, n, E) in Eq. JHJ implies that the luminosity at 
radius r must be interpreted as originating from a neu- 
trino number per proper time crossing a mass shell as 
measured by an observer comoving with the shell. The 
neutrinos carry energies, E, which are also measured in 
the comoving frame. 

We compare time slices in our runs at three crucial 
instances in the postbounce evolution: Bounce (t p b = 
ms), burst (f p (, = 3 ms), and broil (t p i > 100 ms). The 
importance of core bounce, the instance when the in- 
fall is reversed at the center due to the strong repulsive 
forces above nuclear density, needs no further explana- 
tion. The time slice at 3 ms not only captures the launch 
of the electron neutrino burst, but also the interesting 
phase when the dynamical bounce-shock stalls (i.e., the 
postshock velocities become negative). At this time, long 
before any neutrino heating can take place, the infalling 
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material does no longer change the direction of its veloc- 
ity at the shock front. After deceleration at the accre- 
tion shock it continues to drift toward the center of the 
star. Puffed-up by the dissipation of the kinetic energy 
acquired during infall, however, the net volume of accu- 
mulated shock-heated material still increases such that 
the accretion front continues to expand to larger radii 
in a quasi-stationary manner. Only after the accretion 
front has reached a position farther out in the gravita- 
tional well, where less kinetic energy is dissipated in the 
shock, do the temperature difference between the reced- 
ing hot neutrinospheres and the cooler postshock matter 
become favorable for neutrino heating. These conditions 
are not met before a time of 50 ms after bounce. In spher- 
ically symmetric simulations the accretion front reaches 
a maximum radius around 150 km at about 100 ms af- 
ter bounce and recedes slowly thereafter. Therefore, we 
choose a third snapshot in our comparison at this late 
phase where neutrino cooling and heating influence the 
quasi-stationary evolution. 

4.1. Hydrodynamics 

As described above, the dynamical simulations are 
based on the two very different hydrodynamics codes 
agile and PROMETHEUS. In order to disentangle hydro- 
dynamics differences from neutrino transport differences 
in our results it is helpful to perform a comparison of 
an adiabatic collapse where all neutrino interactions are 
suppressed. We found a simple test case that poses sim- 
ilar challenges to the hydrodynamics algorithms as the 
case with full transport. We take the progenitor model 
of run N13 and replace the electron fraction and entropy 
as a function of enclosed mass by the values obtained 
at bounce in the N13 run with full transport. By this 
measure, the Chandrasekhar mass at core bounce, which 
depends on electron fraction and temperature, is imposed 
already at the beginning of the simulation. As expected, 
the adiabatic collapse of this modified progenitor leads to 
bounce and shock formation at a similar mass coordinate 
as in the simulation with full transport, and therefore to 
similar conditions around bounce. 

Figurenshows the situation at 3 ms after bounce. This 
is the critical time when the dynamical shock in the full 
models stalls to turn into an accretion front. We find 
very good agreement in the density and velocity profiles. 
The latter show the same timing and amplitude of re- 
flected sound waves in the ringing neutron star. Also the 
entropy profiles agree well. The entropy profile contains 
information about the evolution of the shock strength be- 
cause the otherwise conserved entropy of a fluid element 
can only change due to the dissipation of kinetic energy 
when matter passes through the shock front. The pro- 
file shows that both the shock in AGILE and the shock 
in Prometheus start with a similar strength (almost 
identical entropy peak at 0.55 M Q ). In the further evo- 
lution, however, the shock in agile decays slightly faster, 
producing lower postshock entropies. The weaker shock 
propagates more slowly with respect to the mass coordi- 
nate so that a small offset of the shock position becomes 
visible at 3 ms after bounce. 

4.2. Newtonian 13 M Q Model 

For the runs that include neutrino transport, we start 
the comparison with an investigation of differences in the 
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Fig. 1. — Snapshots at 3 ms after bounce for the adiabatic 
collapse of the modified N13 progenitor. Data from the AGILE sim- 
ulation arc drawn with thick lines. Data from the PROMETHEUS 
simulation are drawn with thin lines. Panel (a) shows the veloc- 
ity (solid lines) and density (dashed lines) profiles, panel (b) the 
entropy (solid lines) and electron fraction (dashed lines). The hy- 
drodynamics simulations agree well. The shock strength decays 
slightly faster in AGILE than in PROMETHEUS. The weaker shock in 
AGILE tends to propagate more slowly in mass and to produce a 
smaller postshock entropy. 



N13 model, in which only v e and v & are taken into ac- 
count. Figure |21 presents the conditions at bounce. The 
neutrino luminosities are shown in panel (a). We find 
transient differences of order 30% in the electron neu- 
trino luminosities at bounce in the diffusive regime inside 
of the nascent shock. The luminosities are in good agree- 
ment outside of the shock front and the electron neutrino 
rms energies in panel (c) agree well in all domains. The 
electron antincutrino rms energies are very uncertain at 
this time because the antineutrino abundance is negligi- 




M [M ] M [M ] 

Fig. 2. — Snapshots at bounce for model N13. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. Data from the 
vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, panel (d) the 
entropy (solid lines) and electron fraction (dashed lines). The neutrino luminosities are given in panel (a) and the neutrino rms energies 
in panel (c). Solid lines refer to electron neutrinos and dashed lines to electron antineutrinos. The central electron fraction in vertex is 
smaller than in AGILE-BOLTZTRAN and the shock forms at a slightly smaller enclosed mass. The high electron degeneracy during collapse 
suppresses electron antineutrino production so that the corresponding luminosity at bounce is below the threshold of panel (a). 



ble under the electron-degenerate conditions at bounce. 
Important quantities at bounce are the entropy and the 
electron fraction in panel (d) because they determine the 
size of the causally connected homologous core. A shock 
forms when the outgoing pressure wave from the bounce 
at nuclear densities reaches its edge. The differences in 
the electron fraction are of order 3%, the differences in 
the neutrino abundances are even smaller. 



Most of these differences are introduced during the last 
2 ms before bounce. The profiles are in nearly perfect 
agreement before. This becomes evident in panels (c) 
and (d) of Fig. |3 where we plot the entropy and lep- 
ton fractions of the innermost zone as functions of den- 
sity during core collapse. We find that the differences 
in deleptonization appear just before neutrino trapping 
sets in, i.e. when the effective electron capture rates are 
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highest. The entropy evolution shows perfect agreement 
during infall, but after neutrino trapping, an unphysi- 
cal entropy increase in the innermost zones takes place 
in the AGILE-boltztran simulation. Fig. 01 clearly 
demonstrates that this entropy increase only occurs in 
the innermost zone. Probably more significant is the 
small deviation of order 5% in the lepton and electron 
fraction which appears at the same time. It is not con- 
fined to the innermost zone and influences the enclosed 
mass at shock formation. 

The formation of the shock is visible in the velocity 
profiles in panel (b) of Fig. The difference in enclosed 
mass ~ 3% between agile-boltztran and vertex is 
qualitatively consistent with the differences in the elec- 
tron fraction profiles in Fig. [5Ji. The infall velocities 
in the outer core agree well. Panel (b) also shows the 
density profiles of the N13 run at bounce. We conclude 
the discussion of bounce with the observation that there 
are small visible differences between the two N13 simula- 
tions, but none of them is likely to induce large deviations 
in the postbounce evolution. 

Figure shows the mass trajectories for both runs dur- 
ing the first 10 ms after bounce, in panel (a) for AGILE- 
BOLTZTRAN and in panel (b) for vertex. The rising 
dashed line marks the shock position. A first inspec- 
tion reveals a difference of 15% in shock radius at 7 ms 
after bounce. As we will see later in Fig. [SJ this dif- 
ference is transient. It appears after a short dynami- 
cal phase of shock propagation, long before any neutrino 
heating takes place and long before the accretion front 
has reached its maximum radius at ~ 250 km in this 
optimistic model with reduced input physics. 

The discrepancy originates from a different shock 
strength as cause and a different neutrino burst behavior 
as consequence and amplification mechanism. The gra- 
dients of the mass trajectories in Figs. and b indicate 
the velocity of the material in the postshock region (at 
the right hand side of the thick dashed line that repre- 
sents the position of the shock). The bounce-shock is 
dynamical at the beginning and drives matter outward 
with positive postshock velocities. At about 4 ms after 
bounce, the shock stalls because of the nuclear dissocia- 
tion of infalling matter and neutrino losses. It converts 
into an accretion shock, characterized by jump condi- 
tions that connect the high speed/low density infalling 
material to the low speed/high density postshock ma- 
terial. There is an important difference between the dy- 
namical and the accretion shock: The postshock material 
behind a dynamical shock is diluted between the rarefac- 
tion wave and the shock front, while the material behind 
an accretion front continues to be compressed due to the 
accumulated mass. The examination of the mass trajec- 
tories in Fig. [jjib indicates that vertex maintains a dy- 
namical shock for a longer time than agile-boltztran. 
This is consistent with the result of the hydrodynamics 
comparison in Fig. ^ 

This difference in shock propagation is initially not 
very large. It is significantly amplified by the coincidence 
of the electron neutrino burst with the transition from 
a dynamical to an accretion shock. As the shock com- 
presses the infalling lepton-rich material, the fermionic 
electrons have to populate high energy levels and are 
rapidly converted to neutrinos by captures on protons as 
soon as the density is low enough for the produced neu- 



trinos to escape. This neutrino burst can be launched 
by the shock while it is still in its dynamical phase or 
after it has stalled to an accretion front. If the neutrino 
burst is launched during the dynamical phase, infalling 
matter deleptonizes due to the immediate compression 
in the shock front. After that, electron captures are 
reduced quickly because the matter re-expands behind 
the dynamical shock and the density drops again. In 
an accretion shock, the infalling matter experiences the 
same initial deleptonization in the shock front, but con- 
tinues to emit neutrinos due to the compression behind 
the shock. The neutrino burst from an accretion shock 
is therefore more intense. The neutrino emission behind 
the shock removes energy and lepton number and reduces 
the pressure support. It accelerates the attenuation of 
the weaker shock and thereby amplifies the initial differ- 
ence in shock strength. Using the lepton number source 
term, qi, in units of generated leptons per baryon and 
second, we shaded in panels (a) and (b) of Fig. the 
regions where the neutrino emission 4irr 2 p/ms X qi ex- 
ceeds the thresholds of 10 51 erne's" 1 , 2 x 10 51 erne's" 1 , 
and 3 x 10 51 erne's -1 . The comparison of panels (a) and 
(b) illustrates that more neutrinos are lost from the re- 
gion behind the weaker shock in the agile-boltztran 
simulation. The initially stronger shock in the vertex 
simulation suffers less neutrino losses and leads to a more 
rapidly expanding accretion front. 

After this investigation we can easily interpret the time 
slice at 3 ms after bounce presented in Fig. 01 In the 
luminosity, panel (a), we see good agreement during the 
launch of the neutrino burst. The somewhat higher rms 
energies of the burst electron neutrinos, compared to the 
previously emitted ones, are visible in panel (c). The 
entropy profile in panel (d) shows a slightly weaker shock 
in agile-boltztran than in vertex, very similar to the 
differences found for the hydrodynamics comparison in 
Fig. n The electron fraction profiles now show that the 
deleptonization of the postshock region between 0.8 and 
1.0 M Q is more pronounced in agile-boltztran. This 
is consistent with the higher density visible in panel (b). 
The positive velocities behind the shock also demonstrate 
that the vertex shock has still a larger kinetic energy in 
this snapshot. The shock in vertex will stall somewhat 
later than in agile-boltztran. 

The further evolution is best followed in the time- 
dependent diagrams in Fig. [SJ In panel (b) , the luminosi- 
ties and rms energies are shown, sampled at a fixed radius 
of 500 km. The two neutrino signals are qualitatively 
very similar. The neutrino burst in agile-boltztran is 
somewhat broader than in vertex and has a 7% smaller 
peak luminosity. The deviations after the burst are at 
most 15% around 50 ms after bounce. This difference 
is a late consequence of the deleptonization differences 
during the neutrino burst. The material in vertex is 
left with higher electron fraction and higher entropy af- 
ter the neutrino burst. Hence it deleptonizes at a higher 
rate afterwards. The rms energies tend to be lower in 
vertex than in agile-boltztran. Finally, panel (a) 
compares the accretion front trajectories over a longer 
period of time. We find the described differences in the 
early expansion of the accretion front. At 85 ms after 
bounce, however, the trajectories cross and a larger max- 
imum radius is reached in the agile-boltztran simula- 
tion. The maximum radius of the accretion front differs 
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Fig. 3. — Panel (a) shows the first 10 ms after bounce for model N13 in the simulation with AGILE-BOLTZTRAN. Panel (b) shows the 
same time period in the simulation with vertex. The thin lines represent trajectories of fluid elements spaced with an interval of 0.02 
Mq. The dashed line marks the shock position as a function of time. The gradients of the mass trajectories at the right hand side of the 
dashed line indicate the postshock velocities. At ~ 4 ms after bounce the shock has stalled. It turns into an expanding accretion front 
with negative postshock velocities. Areas with strong neutrino emission are shaded in three levels corresponding to values of one, two, 
and three times 10 51 neutrinos per centimeter and second (i.e. for Airr 2 p/rriB X qe, where qi is the lepton number source term in units of 
leptons per baryon and second). The coincidence of the launch of the neutrino burst with the transition from a dynamical to an accretion 
shock in this model leads to an amplification of the small hydrodynamics differences found above. The region behind the weaker shock 
in AGILE-BOLTZTRAN experiences more compression. It therefore deleptonizes more rapidly and the weak shock loses even more pressure 
support than the stronger shock in the vertex simulation. The accretion front therefore expands more slowly in the AGILE-BOLTZTRAN 
simulation. Panel (c) compares the entropy of the innermost zone as a function of density in AGILE-BOLTZTRAN (thick line) and in VERTEX 
(thin line). The agreement before trapping is close to perfect, dynamically insignificant differences appear at larger densities. Panel (d) 
shows an analogous comparison for the electron fraction (solid lines) and lepton fraction (dashed lines) in the innermost zone. 
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Fig. 4. — Snapshots at 3 ms after bounce for model N13. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. Data 
from the vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, panel 
(d) the entropy (solid lines) and electron fraction (dashed lines). Panels (a) and (c) show the neutrino luminosities and rms energies 
respectively. Solid lines refer to electron neutrinos and dashed lines to electron antineutrinos. The snapshot reveals similar small differences 
as the hydrodynamics comparison in Fig. 
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Fig. 5. — The shock position as a function of time for model N13 is shown in panel (a). The shock in vertex (thin line) propagates 
initially faster and nicely converges after its maximum expansion to the position of the shock in AGILE-boltztran (thick line). The neutrino 
luminosities and rms energies for model N13 are presented as functions of time in panel (b). The values are sampled at a radius of 500 
km in the comoving frame. The solid lines belong to electron neutrinos and the dashed lines to electron antineutrinos. The line width 
distinguishes between the results from AGILE-BOLTZTRAN and vertex in the same way as in panel (a). The luminosity peaks are nearly 
identical, the rms energies have the tendency to be larger in AGILE-BOLTZTRAN. 



by ~ 8% between the two simulations. This difference 
is due to the higher preshock entropies in the AGILE- 
BOLTZTRAN simulation visible in Fig. EJi. The differ- 
ence in the entropies of the infalling material stems from 
the interface between the silicon layer and the material 
in nuclear statistical equilibrium (NSE). The burning in 
AGILE-BOLTZTRAN cannot be switched off completely for 
this comparison, because conversions between non-NSE 
and NSE are unavoidable when the adaptive grid moves 
its zones relative to the mass coordinate. This produces 
a local entropy peak at the composition interface (see 
Fig. 0Jl at an enclosed mass of 1.2 M ) which crosses 
the accretion front at 85 ms after bounce. The higher en- 
tropy leads to a temporarily lowered accretion rate which 
allows the accretion front in the AGILE-BOLTZTRAN sim- 
ulation to propagate to a larger radius than in the VER- 
TEX simulation where no conversions between non-NSE 
and NSE occur during this simulation with reduced input 
physics. After the initial expansion phase, where matter 
piles up on the neutron star hydrostatically, the pressure 
support in the cooling region starts to diminish rapidly 
and the matter in t he heat i ng reg i on is pulled inward 
from b elow (see e.g. IJanka I 1120011) ; ILiebendorfer et al~l 
1120011) ). Simultaneously, the mass shell in the preshock 
region containing the described entropy differences has 
fallen through the shock and the entropy becomes more 
similar again. Therefore the trajectories of the accretion 
front converge and agree well during the shock recession 
phase. 

We finish the comparison of the N13 model with a 
closer look at a time slice at 150 ms after bounce (Fig. 



|SJ), which is a snapshot during this quasi-stationary ac- 
cretion phase. Panel (a) demonstrates excellent agree- 
ment of the luminosities in all regions of the computa- 
tional domain. The rms energies in panel (c) show small 
differences, especially outside of 100 km radius. The ve- 
locity profiles in panel (b) agree well if one disregards the 
different shock positions explained above. The higher ac- 
cretion rates in vertex of material with lower entropy 
lead to a higher density in the cooling region. This is vis- 
ible in the density profiles in panel (b) and the entropy 
profiles in panel (d). The reaction time scale is compa- 
rable to the dynamical time scale inside the gain radius 
at 115 — 120 km. The infalling fluid element therefore 
is close to weak equilibrium in the given neutrino back- 
ground. Since the neutrino luminosities are very similar 
in the two simulations, the larger density of the vertex 
run requires a lower equilibrium electron fraction in the 
cooling region inside the gain radius. The corresponding 
differences in the pressure profiles imply less support for 
the shock and are consistent with a smaller radius of the 
accretion front. 

4.3. General Relativistic 15 Mq Model 

For the analysis of the G15 simulations we start again 
with the description of the conditions at bounce. The 
evolution of the entropy and lepton fraction in the in- 
nermost zone is shown in panels (c) and (d) of Fig. [7| 
Compared to model N13 with Newtonian gravity, the 
central entropy in the general relativistic G15 model is 
15% higher and the lepton and electron fractions are 5% 
lower. The deviations between the vertex and agile- 
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Fig. 6. — Snapshots at 150 ms after bounce for model N13. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. Data 
from the vertex simulation arc drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, panel 
(d) the entropy (solid lines) and electron fraction (dashed lines). The neutrino luminosities and rms energies are shown in panels (a) and 
(c) respectively. Solid lines refer to electron neutrinos and dashed lines to electron antineutrinos. This stationary-state situation is typical 
of the neutrino heating phase at later time after bounce. The agreement is satisfactorily close in most quantities. 



BOLTZTRAN simulations are of order 3% with the excep- 
tion of an entropy increase in the innermost zone in the 
AGILE-BOLTZTRAN simulation. As shown in Fig. IBJi, 
excellent agreement is found in all other regions of the 
model up to the burning front, where different approxi- 
mations in tracking the composition and nuclear burning 
explain the larger differences. In contrast to the entropy 
difference in the innermost zone, the small differences 
in the electron fraction apply to the whole high-density 



region enclosed by the shock. 

The luminosity profiles at bounce are displayed in Fig. 
IHt- Outside of the shock front, they have been set dur- 
ing collapse and agree well. Also the luminosities in the 
diffusive regime reveal no mentionable differences. As a 
consequence of the differences of the electron fraction vis- 
ible in panel (d), the velocity and density profiles in panel 
(b) show the formation of the shock front in vertex at a 
slightly deeper point than in AGILE-BOLTZTRAN. In con- 
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Fig. 7. — Panel (a) shows the first 10 ms after bounce for model G15 in the simulation with AGILE-boltztran. Panel (b) shows the 
same time period in the simulation with vertex. The thin lines represent trajectories of fluid elements spaced with an interval of 0.02 
Mq. Areas with strong neutrino emission are shaded in three levels corresponding to values of one, two, and three times 10 51 neutrinos 
per centimeter and second (i.e. for Airr 2 p/mg X qi, where qi is the lepton number source term in units of leptons per baryon and second). 
Both codes obtain an extended region of strong neutrino emission behind the shock, which turns into an accretion front (dashed line) at 
3 — 4 ms after bounce. Panel (c) compares the entropy of the innermost zone as a function of density in AGILE-BOLTZTRAN (thick line) and 
in vertex (thin line). The agreement before trapping is close to perfect, dynamically insignificant differences appear at larger densities. 
Panel (d) shows a comparison of the electron fraction (solid lines) and lepton fraction (dashed lines) in the innermost zone. The deviations 
between the two codes are of order 3%. 



sideration of this displacement in the shock position, also 
the rms neutrino energies in panel (c) are in satisfactory 
agreement. 

The early postbounce evolution of the G15 model 
is less sensitive to small differences than the pre- 
viously discussed N13 simulation. The main rea- 
son is the weaker bounce shock. General relativis- 



tic effects during core collapse shift the sonic point 
to a 20% smaller enclose d mass and lead to a 
less energetic bounce shock llLiebendorfer et al. I I200H 
iBruenn. DeNisco. fc Mezzacanna I 120011) which has to 
dissociate more infalling material. The inclusion of /i- 
and r-neutrinos in the G15 runs causes additional en- 
ergy drain from the region behind the shock. After the 
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Fig. 8. — Snapshots at bounce for model G15. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. Data from 
the vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, panel (d) 
the entropy (solid lines) and electron fraction (dashed lines). The neutrino luminosities and rms energies are given in panels (a) and (c), 
respectively. Solid lines correspond to electron neutrinos, dashed lines to electron antineutrinos, and dash-dotted lines to fi- or r-neutrinos 
(or their antiparticles). In the region enclosed by the shock the luminosities in AGILE-BOLTZTRAN are smaller than in vertex, consistent 
with the larger central electron fraction and a slightly larger enclosed mass at shock formation. 



shock has stalled within 1 ms in both G15 simulations, 
the electron neutrino burst is launched during the ac- 
cretion shock phase. As expected from the discussion of 
the neutrino emission from the shock in the N13 model, 
panels (a) and (b) in Fig. reveal a well extended re- 
gion of high neutrino emission behind the shock. Since 
the postshock matter develops negative velocities a few 
milliseconds after bounce, the further evolution is deter- 



mined by the continued accumulation of accreted mat- 
ter, more and more effectively cooled as the accretion 
front reaches layers with lower matter densities and neu- 
trino opacities. This quasi-stationary evolution is less 
sensitive to differences in the numerics or input physics 
than the dynamical shock propagation in the more opti- 
mistic N13 simulation. The feedback between neutrino 
transport and hydrodynamics amplifies differences in the 
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propagation of dynamical shocks, because the material 
behind the weaker shock emits more neutrinos so that 
the shock loses even more pressure support. In case of 
a quasi-stationary accretion front, differences are likely 
to be reduced because a larger accretion rate produces 
larger neutrino losses. 

However, differences can still be observed: The AGILE- 
BOLTZTRAN shock is stronger at formation. As in the 
N13 simulations, this is partly due to the higher electron 
fraction in the homologous core of the AGILE-BOLTZTRAN 
simulation. A difference in the entropy profiles is left 
behind when the shock passes an enclosed mass of 0.75 
M© in Fig. EJi. This points to a larger deviation between 
the initial shock strengths than in the N13 simulations 
so that the AGILE-boltztran shock expands initially 
faster in Fig. [7^,. The differences between the electron 
fraction profiles in Fig. however, suggest that later 
on a faster deleptonization in AGILE-boltztran damps 
the expansion of the accretion front and as a consequence 
the shock positions in both simulations converge again. 
This connection between a somewhat enhanced neutrino 
loss and a deceleration of the expansion of the accretion 
front is supported by the shaded areas in panels (a) and 
(b) of Fig. [7| which highlight differences in the regions 
of strong neutrino emission in both simulations. 

In Fig. Eb, a lower density at 0.75 M© is caused by 
the higher entropy in the AGILE-boltztran run. The 
velocity profiles at 3 ms after bounce are in very good 
agreement. A slightly higher infall velocity in front of 
the shock leads to a positive entropy gradient between 
0.8 M© and ~ 1 M© in vertex. Also the neutrino lumi- 
nosities at 3 ms after bounce in Fig. [S^ do not reveal new 
features, except for the presence of fj,- and r-neutrinos 
that had not been included in the N13 run. Because of 
the absence of charged-current reactions of these neu- 
trinos, they decouple at higher densities and reach ap- 
preciable luminosities earlier in the evolution, before the 
appearance of electron antineutrinos, which are initially 
suppressed by the high electron degeneracy. The appar- 
ent difference in the electron antineutrino luminosity be- 
tween AGILE-BOLTZTRAN and vertex is due to a small 
time lag in displaying this rapidly rising quantity. 

The further evolution is resumed in Fig. 1101 Panel 
(a) compares the position of the accretion front as a 
function of time for the G15 models. During the expan- 
sion of the accretion front, we find very good agreement 
and the maximum radius is nearly identical. Afterwards, 
the accretion front retreats somewhat more slowly in the 
AGILE-boltztran simulation than in the vertex simu- 
lation. In the latter, the retraction transiently stagnates 
between 150 and 180 ms after bounce. Such features are 
caused by the steep density drop at the infalling inter- 
faces between layers of different composition outside of 
the iron core. The transition to the oxygen-rich silicon 
layer passes the shock at about 165 ms after bounce. 
AGILE-BOLTZTRAN tracks the structure of the outer lay- 
ers less accurately because of the artificial diffusion intro- 
duced by the adaptive grid. Discrete transitions between 
layers are therefore washed out to some extent such that 
their impact on the trajectory of the accretion front is 
less pronounced, although still qualitatively visible. The 
luminosity peaks during the electron neutrino burst de- 
viate only by 3% in the G15 simulations, the average 
peak value is 3.8 x 10 53 erg/s with a half- width of 6 ms. 



The further time evolution of the luminosities and rms 
energies in panel (b) reveals ~ 20% larger values in the 
vertex run. These are at least in part a consequence 
of the increased accretion rate in the vertex simulation 
during the retraction phase of the accretion front. We 
will make this argument more precise in the following 
discussion of late time slices. 

We present two time slices for the long-term evolution 
of the G15 simulation. The first time slice is at 100 ms 
after bounce when the neutrino heating is most efficient. 
The second time slice at 250 ms marks the end of the 
time period covered by the simulations. The time slice 
at 100 ms is given in Fig. Panel (b) shows the shock 
in vertex at a slightly smaller radius and the preshock 
infall velocities to be somewhat higher than in AGILE- 
BOLTZTRAN. The luminosities in panel (a) are up to 20% 
larger in vertex. The luminosity discontinuity across 
the shock front caused by the Doppler frequency shift 
and angular aberration for an observer in the comoving 
frame is also larger in vertex because of the larger lumi- 
nosity and the larger velocity jump visible in panel (b). 
Otherwise the relative differences between the two runs 
are just inverse to the situation we analyzed in Fig. El in 
case of the N13 simulation. Now AGILE-boltztran has 
a higher density in the shocked material in panel (b) and 
a correspondingly lower entropy and electron fraction in 
panel (d). Also the neutrino rms energies in panel (c) 
are now lower than in the vertex simulation. 

We also include the latest time slice at 250 ms after 
bounce in Fig. El The panels show the same qualita- 
tive features we have discussed in the context of Fig. 1111 
but at this late time to a much larger extent. The pro- 
toneutron star in the vertex simulation is more com- 
pact, causing the accretion front to reside at a smaller 
radius and the luminosities of all neutrino flavors to be 
larger and to have harder spectra than in the AGILE- 
BOLTZTRAN simulation. Due to the smaller radius of the 
accretion front and higher infall velocities ahead of it, the 
postshock entropy is significantly higher in the vertex 
run. 

A more compact neutron star with a larger mass and 
a stronger gravitational potential (which would lead to 
higher preshock velocities) can be a consequence of a 
higher mass accretion rate. Indeed we observe differences 
in the mass flux to the shock, which can be caused by the 
different quality to follow structures in the outer stellar 
layers in both simulations. An infalling density feature 
in vertex (which might evolve differently due to the dif- 
ferent treatment of "burning" and thus different entropy 
and pressure, or may be smoothed by the diffusivity of 
the adaptive grid in AGILE-boltztran) transiently en- 
hances or reduces the accretion rate. When, for exam- 
ple, the transition to the oxygen-rich silicon layer falls in 
shortly before 200 ms after bounce, the accretion rate de- 
creases sharply and the retraction of the accretion front 
in the vertex simulation stagnates. This is consistent 
with the luminosity reduction during this phase. But 
the artificial diffusion in the outer layers of the AGILE- 
BOLTZTRAN simulation can explain only transient differ- 
ences of the mass accretion rate and of the total accreted 
mass, because an associated redistribution of matter and 
a modification of the preshock structure is limited to a 
certain radial domain. It should, however, not produce 
persistent differences in the density distribution behind 
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Fig. 9. — Snapshots at 3 ms after bounce for model G15. Data from the AGILE-boltztran simulation are drawn with thick lines. Data 
from the vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, panel 
(d) the entropy (solid lines) and electron fraction (dashed lines). The neutrino luminosities and rms energies are given in panels (a) and (c), 
respectively. Solid lines correspond to electron neutrinos, dashed lines to electron antineutrinos, and dash-dotted lines to fi- or r-neutrinos 
(or their antiparticles). Variations in the entropy profiles reflect the differences in the shock strength at bounce. The deleptonization by 
the neutrino burst occurs an instant earlier in the AGILE-boltztran simulation. 



the shock, which in fact can be seen in Figs. 1111 and 1121 
(cf. the density profiles in panels (b)). A closer inspec- 
tion of our models in fact reveals that the mass accretion 
rate outside of the shock and the baryonic mass accumu- 
lated in the neutron star show temporary differences only 
between ~30 ms and ^200 ms, but become very similar 
again toward the end of our simulations. 

The systematically evolving and growing differ- 
ence during the long-term evolution must therefore 



be caused by another effect. The combination of 
higher luminosities, higher rms energies and higher en- 
tropies at the neutrinosphere reminds us of the dif- 
ferences found between Newtonian and general rel- 
ativistic simulations, where they are due to differ- 
ences between the Newtonian and relativistic gravita- 
tional potential ( |Bruenn 1 DeNisco. fc Mezzacappa 120011: 
iLiebendorfer et aLTl^OOlf) . Both of our numerical meth- 
ods are designed to accurately describe the hydrostatic 
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Fig. 10. — Panel (a) shows the position of the accretion front as a function of time for model G15. The two simulations predict a 
very similar early expansion of the accretion front up to a maximum radius of 150 km. Then, the accretion front in vertex (thin line) 
retreats faster than in AGILE-BOLTZTRAN (thick line) in response to the slightly stronger contraction of the nascent neutron star due to the 
approximate treatment of general relativity in the vertex simulation. The hump in the vertex shock position at ~ 180 ms after bounce 
corresponds to an entropy and density discontinuity at the bottom of the oxygen-rich silicon shell. Because the vertex simulation resolves 
the structure of the outer layers of the progenitor star more accurately than the diffusive adaptive grid in AGILE-BOLTZTRAN, this feature is 
less pronounced in the latter simulation. The neutrino luminosities and rms energies are presented as functions of time in panel (b). The 
values are sampled at a radius of 500 km in the comoving frame. The solid lines correspond to electron neutrinos, the dashed lines to electron 
antineutrinos, and the dash-dotted lines to fi- or r-neutrinos (or their antiparticles). The line width distinguishes between the results from 
AGILE-BOLTZTRAN (thick lines) and vertex (thin lines). The differences in the neutrino results are mainly — but not exclusively — indirect 
consequences (due to the more compact neutron star) of the approximate treatment of general relativity in the vertex simulation. 



structure of the protoneutron star according to the so- 
lution of the Tolman-Oppenheimer-Volkoff (TOV) equa- 
tion. But in the relativistic case differences in the po- 
tential can not only be caused by a difference of the en- 
closed mass at a given radius. In contrast to the Newto- 
nian case, the relativistic potential depends highly non- 
linearly on the structure of the configuration through its 
dependence on the mass distribution, pressure, and en- 
ergy density. The gravitational potential is therefore sen- 
sitive to differences in the early post-bounce dynamics of 
the propagating shock (e.g., due to the different initial 
shock strength) and to pressure and entropy differences 
created at later times, e.g., associated with the transient 
differences of the mass accretion rate or due to the higher 
infall velocities ahead of the shock in the vertex run. 
The overestimation of the velocities in the collapse layer 
is also a consequence of the approximation of general rel- 
ativity in case of vertex. The latter code uses only a 
gravitational potential that is corrected for general rela- 
tivistic effects, but ignores relativistic kinematics. When 
the preshock values of the infall velocities reach 10-15% 
of the speed of light, the velocities computed by vertex 
are overestimated in comparison to the relativistic ve- 
locities calculated by agile-boltztran. This again has 
an influence on the long-term post-bounce evolution and 
thus feeds back into the core structure and causes a non- 
linear response of the relativistic potential. While we find 



that the maximum densities at bounce and the following 
relaxation to a static situation are in good agreement be- 
tween the vertex and AGILE-BOLTZTRAN runs, we sub- 
sequently observe clear deviations of the central density 
and of the density profile which gradually evolve and 
grow at later times. The vertex simulation develops a 
higher central density and a steeper density gradient out- 
side of the high-density core, and thus a more compact 
neutron star with a higher relativistic potential. This 
is confirmed by Fig. 1131 which shows the profiles of the 

metric coefficients a = g\[ 2 and T — g a a dr/da in Eq. 

at 250 ms after bounce. The smaller deviations from 
unity of the metric coefficients in the AGILE-BOLTZTRAN 
run are consistent with the less compact structure of the 
protoneutron star in this simulation. The profiles also 
show that the metric coefficients are nearly unity outside 
of the accretion front. The larger preshock velocities for 
the VERTEX run in Figs. Illb and 1 12b are therefore not a 
consequence of the fact that the lapse function is not in- 
cluded in the velocity plotted for this simulation. They 
are more likely caused by the disregard of kinematical 
effects and a stronger gravitational potential in the ap- 
proximation of general relativity used in vertex. 

4.4. Discussion 

This work extends the testing of our codes 
beyond the independently performed calcula- 
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Fig. 11. — Snapshots at 100 ms after bounce for model G15. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. 
Data from the vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, 
panel (d) the entropy (solid lines) and electron fraction (dashed line). The neutrino luminosities and rms energies are displayed in panels 
(a) and (c), respectively. Solid lines correspond to electron neutrinos, dashed lines to electron antineutrinos, and dash-dotted lines to fx- 
or t- neutrinos (or their antiparticles). In this time slice we find a smaller shock position, somewhat faster infall ahead of the shock, and 
higher postshock infall velocities in vertex compared to the results of AGILE-BOLTZTRAN. Consistent with the more compact structure of 
the protoneutron star, the entropies, neutrino luminosities, and neutrino rms energies in vertex are larger than in AGILE-BOLTZTRAN. 



tions of idealize d probl e ms w hic h have analyt- 
ical s olutions iMes ser I 12 0001 IRampp fc Janka I 
1 20021 lUebendorfer. Rosswog. k. Thielemann I I2002t 
iLiebendorfer et al. I l2004jh Here we directly compare 
the results of the codes in the application they were ac- 
tually developed for. Our aim was to assess quantitative 
differences in complex supernova simulations to reduce 
the probability of qualitative differences in future appli- 



cations. We also intended to create points of reference 
for future testing of codes that handle the challenges 
of supernova physics, and to lay the foundations for 
performing such tests in a more realistic way than by 
means of comparison with analytic solutions of idealized 
test problems. 

We have encountered two fundamental difficulties in 
our comparison. The first is the fact that the two codes 
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Fig. 12. — Snapshots at 250 ms after bounce for model G15. Data from the AGILE-BOLTZTRAN simulation are drawn with thick lines. 
Data from the vertex simulation are drawn with thin lines. Panel (b) shows the velocity (solid lines) and density (dashed lines) profiles, 
panel (d) the entropy and electron abundances. The neutrino luminosities and rms energies are displayed in panels (a) and (c). Solid lines 
correspond to electron neutrinos, dashed lines to electron antineutrinos, and dash-dotted lines to /i- or r-ncutrinos (or their antiparticles) . 
The protoncutron star in the vertex simulation is more compact than in the AGILE-BOLTZTRAN simulation. 



employ different methods, use different basic quantities, 
and are differently structured. The comparison of the 
results of the two approaches is straightforward, but it 
is by far more challenging (and sometimes impossible) 
to track differences back to their origin if the compared 
quantities are not calculated in a similar way. The sec- 
ond difficulty is more related to supernova physics than 
to the methodology. In the comparison of our results we 
have very often encountered the situation that all quanti- 



ties within either one simulation are perfectly consistent, 
but still not the same as in the other simulation. Further 
investigations of the differences revealed small initial dif- 
ferences in several tightly coupled quantities that grow 
with ongoing evolution. Because of the strong feedbacks 
in the supernova problem it was often almost impossible 
to separate cause and consequences of the deviations. 
This strong coupling between quantities indicates that 
the problem is governed by many equilibria. Sometimes 
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g a a /2 dr/da 



the results converge again and find back to a similar evo- 
lution as soon as the equilibrium is achieved. This is for 
example the case in model N13 when stationary-state 
accretion is established after the very dynamical early 
postbounce phase. 

When we consider just the overall dynamical evolution 
of the spherically symmetric postbounce phases in Figs. 
and 1101 we find agreement in all qualitative features 
of the history. But we also find some significant quanti- 
tative differences. For example, the early shock propa- 
gation and the luminosities during the first 100 ms after 
bounce are different in model N13. We attribute an im- 
portant part of these differences to the hydrodynamical 
shock propagation that appears to maintain a stronger 
shock in vertex than in AGILE-BOLTZTRAN. The close 
coincidence of the neutrino burst with the transition from 
a dynamical to an accretion shock in this model leads to 
an amplification of existing differences. The region be- 
hind the weaker shock in AGILE-BOLTZTRAN deleptonizes 
to a larger extent than behind the slightly stronger shock 
in vertex. While the deleptonization burst is more 
extended in AGILE-BOLTZTRAN, the vertex core could 
maintain higher luminosities later on. This may be the 
reason for a somewhat more optimistic shock propaga- 
tion in the early phase of the vertex simulation of model 
N13. 

No such amplification effect takes place in the G15 
model where the shocks in both runs have already made 
the transition to an accretion front when they reach den- 
sities from where neutrinos begin to escape. Despite 
of some differences in the time-dependence of the shock 



strength, the evolution of these two more realistic runs 
with "standard" input physics agrees nicely during the 
early post-bounce phase. We find good agreement of 
the neutrino quantities in the diffusive inner core of the 
protoneutron star and the timing and peak height of 
the neutrino burst. The approximations of general rel- 
ativistic effects in vertex yield accurate results until 
bounce and do not introduce larger uncertainties with 
respect to the general relativistic approach than other 
acceptable approximations. Differences appear only in 
the later evolution when the outer layers of the progen- 
itor fall into the stalled accretion front. Some of these 
differences are caused by a different description of nu- 
clear burning in both codes and a different capability to 
track the composition interfaces in the outer layers of 
the collapsing core. The main reason for the differences, 
however, is the approximate treatment of general rela- 
tivity in the vertex simulation. The basically Newto- 
nian hydrodynamics code employs a relativistically mod- 
ified gravitational potential which in principle allows one 
to accurately describe hyd rostatic configurations a ccord- 
ing to the TOV equation ijRampp fc Janka 1 120021 . But 
the code disregards the effects of relativistic kinemat- 
ics which causes an overestimation of the infall velocity 
ahead of the shock. The corresponding differences in 
the dynamical evolution feed back into the gravitational 
potential in a nonlinear way. This leads to a slightly 
more compact neutron star, a somewhat smaller radius of 
the accretion front, and a faster infall of matter between 
shock and neutron star, producing up to ~ 20% higher 
accretion luminosities and rms energies of neutrinos and 
antineutrinos of all flavors. While most of these discrep- 
ancies in the neutrino quantities are a consequence of the 
different structure of the accretion layer in the vertex 
and AGILE-boltztran simulations, a smaller contribu- 
tion may also be ascribed to the fact that vertex takes 
into account general relativistic redshift, but ignores the 
metric effects in the radial coordinate. 

The overall evolution, however, is consistent between 
both runs also in case of the relativistic G15 model. Even 
more, both computations produce not only a remarkable 
qualitative similarity of the behavior during all phases 
but also show nice agreement in most features of the 
radial profiles of the important quantities. This result 
may be especially useful for multi-dimensional simula- 
tions, where essential general relativistic effects should 
not be ignored, but a full relativistic treatment might 
not have highest priority. 

Both methods have their vulnerabilities and some of 
them have led to lively discussions in the past. AGILE- 
BOLTZTRAN was criticized because of the rigorous ap- 
proach and the generous consumption of computer mem- 
ory and CPU time. The resolution of the neutrino phase 
space was considered to be at the lower justifiable limit 
in earlier simulations. And indeed, the certainty that 
the evolution follows basic physical principles indepen- 
dently of the resolution must be earned by very specific 
twists and wrinkles in the finite difference representation 
l|Liebendorfer et al. ll2004|) . vertex has raised concerns 
with regard to consistency as it uses disjunct gridding for 
radiation transport and hydrodynamics and applies re- 
gridding procedures after bounce. The separate solution 
of the transport equations for neutrino number and neu- 
trino energy only adds to the complexity. However, we 
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did not find discrepancies in our comparison that would 
support any of these concerns to a degree that would 
question the reliability of the qualitative results of our 
explosion-free supernova simulations. 

The different strengths of the codes are more visible in 
the quantitative details of the calculations, vertex pro- 
duces great angular resolution in the flux factors even 
far from the neutrinospheres and keeps properly track 
of the sharp discontinuities in the composition of the 
outer layers. Its e xtendibility to two-dim ensional sim- 
ulations is built-in (Ra mpp fe Janka Il2002j) and the gen- 
eral relativistic approximation can be expected to pro- 
duce good results in multi-dimensional situations as well 
ijBuras et al. 1 12003). However, the solution of a model 
Boltzmann equation is much more involved in more than 
one dimension. Additional approximations introduced 
by spherical averaging or ray-by-ray techniques cannot 
be tested in a comparison between spherically symmetric 
models. If these approximations are good, the variable 
Eddington factor approach is a very efficient technique 
for multidimensional simulations. 

AGILE-BOLTZTRAN demonstrates that the solution of 
only one transport equation for the neutrino distribution 
function can provide accurate radiation transport solu- 
tions in the diffusion limit and semi-transparent regime. 
Number and energy conservation are reasonably well ful- 
filled and the description of hydrodynamics and radiation 
transport add up to one consistent general relativistic fi- 
nite difference represent ation of radiation hydrodyn amics 
in spherical symmetry ( Lie bendorfer et al. l l2004'l. The 
approach is in principle extendible to two and three di- 
mensions and allows for adaptive zoning because it is 
based on the local description of the transport equation 
in confined fluid elem ents. But consistency is easily lost 
in higher dimensions ijCardall fe M ezzacappa 2003) an d 
computer performance may become prohibitive in an im- 
plicit multidimensional discrete ordinates approach. 

5. CONCLUSIONS 

We have compared two different approaches to imple- 
ment Boltzmann neutrino transport in spherically sym- 
metric radiation hydrodynamics simulations of stellar 
core collapse and postbounce evolution. We performed 
calculations for two different progenitor stars , the 13 
M© progenitor of lNomo to fc Hashimoto I ill 9881) a nd the 
15 M© progenitor of [Wooslev fc Weaver I 1)1995). We 
present one Newtonian calculation (N13) with the mini- 
mum input physics that leads to a plausible scenario after 
bounce and a second relativistic calculation (G15) with 
the "standard" physics used in many recent supernova 
simulations. We find similar agreement in both cases. 
The reduced complexity of the input physics in the N13 
model helps to isolate differences in the implementation 
of the hydrodynamics and the neutrino transport. We 
could improve the agreement by upgrading the first order 
donor-cell advection scheme in the implicit hydrodynam- 
ics code AGILE to a second order TVD advection scheme. 
The version with first order advection led to a more pes- 
simistic shock propagation during the first 10 ms after 
bounce. It did, however, reveal an interesting relation- 
ship between the transition of the propagating disconti- 
nuity from a dynamical shock to an accretion front and 
the almost coincident launch of the neutrino burst. 

A neutrino burst radiated from an accretion front 



maintains a high luminosity for a longer time than a 
neutrino burst produced by a dynamical shock, because 
an accretion front compresses matter at steady-state like 
conditions whereas the layer behind a dynamical shock 
gets diluted quickly so that electron captures diminish on 
a short timescale. Therefore less lepton number is lost in 
neutrinos from a dynamical shock which rapidly crosses 
the neutrinospheres, but neutrinos extract more leptons 
from the compressed matter behind the accretion front 
once the shock has stalled (i.e., the postshock velocity has 
become negative). This effect, however, turned out to 
produce transient differences only for a few milliseconds 
in our simulations, and convergence of the shock trajec- 
tories was found again after the shocks in both runs had 
transformed to accretion fronts. While the optimistic 
N13 model with only one neutrino flavor (z/ e and D e ) rep- 
resents a case where the neutrinospheres are crossed by 
a dynamical shock, the relativistic model G15 serves as 
an example where the shock forms at a smaller enclosed 
mass due to the deeper general relativistic potential and 
where additional losses occur by the emission of [i- and 
r-neutrinos from deeper layers. In this case the shock 
turns into an accretion front before or at the time the 
neutrino burst is launched. 

The overall evolution of both models is in good agree- 
ment when simulations with the two codes are compared. 
Differences in details were found, e.g., a slightly differ- 
ent shock propagation in the early hydrodynamical phase 
and more smearing of the composition interfaces in the 
outer progenitor layers by artificial diffusion in the case 
of AGILE. The luminosities in vertex tend to be slightly 
higher than in AGILE-BOLTZTRAN and the rms energies 
a little lower in the N13 model. The approximation of 
general relativistic effects by a modified gravitational po- 
tential in otherwise Newtonian hydrodynamics in ver- 
tex is very accurate up to bounce. In comparison with 
the general relativistic simulation of AGILE-boltztran, 
however, a somewhat deeper potential associated with 
higher accretion rates develops during the long-term 
postbounce evolution. The consequence are larger neu- 
trino luminosities and rms energies. But in general, good 
qualitative and satisfactory quantitative agreement of all 
important temporal and radial features was found also 
in the relativistic model. Major differences can result 
from implementation-specific rather than from method- 
specific details, e.g. from the former use of a low-order 
advection scheme in AGILE-boltztran or from the spe- 
cific choice of the finite differencing in both codes. 

We come to the conclusion that both methods work 
satisfactorily well in this application and give compa- 
rable results. We determined similar computational 
needs for our not thoroughly optimized codes. Standard 
runs with AGILE-boltztran tend to consume slightly 
less computer time. But standard runs with vertex 
have been performed with better energy resolution and 
the angular resolution that can be achieved at larger 
radii is out of reach for Sat methods. Hence, a de- 
tailed comparison of CPU time requirements is not re- 
ally meaningful. Moreover, fast er methods may have 
been developed in the meantime <|Burrows et al. 1120001 
iThomoson. Burrows, fc Pinto 1120031) . Rather than argu- 
ing about the "best" method for a certain application, 
we recommend to pursue a variety of feasible numerical 
approaches for future astrophysical simulations, opening 
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up the possibility of independent mutual validation of the 
results. We hope that our comparison provides a useful 
step towards quantitative modeling of a very complex 
astrophysical problem. 
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